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Abstract 


The stability problem of the GEOS satellite has been solved. A com- 
puter simulation indicates lack of stability, a fact that can be attributed 
to the lack of bending stiffness of the cables. Whereas small cable bending 
stiffness can render the system stable, the first natural frequency of 
oscillation of the spacecraft is likely to be very low, so that the cables 
can represent a potential problem area. 



Introduction 


The GEOS satellite (the "simple model") consists of a rigid core, one 
pair of radial booms, one pair of cables with tip masses, and two pairs of 
axial booms, as shown in Fig. 1. The latter two pairs of booms are not 
strictly axial, as they are inclined with respect to the equatorial plane 
at angles other than 90°. The satellite spins freely in space with constant 
angular velocity n. The interest lies in the stability of motion when the 
spacecraft is perturbed slightly from the uniform spin equilibrium state. 

The stability of force-free satellites with flexible appendages, such 
as that considered here, has been investigated on several previous occasions 
(Refs. 1, 2, 3). Such systems are described by both ordinary and partial 
differential equations and are referred to as hybrid . The formulation pre- 
sented in Refs. 1, 2, 3 is perfectly valid for the GEOS satellite. Hence, 
we shall dispense with the details and only outline the method of approach. 

It was shown in Refs. 1, 2, 3 that, under certain circumstances, the 
Liapunov direct method with the Hamiltonian as a Liapunov functional can 
be used to test the stability of hybrid dynamical systems. The main prob- 
lem is how to treat continuous elastic members. The Liapunov direct method 
has been used widely in conjunction with discrete systems. To test the 
stability of an equilibrium point, the testing function must satisfy one 
of several stability or instability theorems. If it does, then the testing 
function is said to be a Liapunov function and appropriate stability, or 
instability, conclusions can be drawn. If it does not, the analysis is 
inconclusive. The stability analysis consists of testing the sign proper- 
ties of the testing function. The problem in applying the Liapunov direct 
method to hybrid systems lies in the difficulty of testing the sign proper- 
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ties of the testing functional (as opposed to the testing function). The 
author of this report has developed and used three different approaches to 
treat the problem of hybrid systems, namely, (1) the method of testing density 
functions, (2) the method of integral coordinates, and (3) the assumed 
modes method. The method of testing density functions works directly with 
the hybrid dynamical system but is quite often unduly restrictive. On the 
other hand, the remaining two methods are based on discretization schemes, 
which implies that the testing functional is replaced by a testing function. 

In particular, the method of integral coordinates involves the definition 
of new generalized coordinates representing certain integrals appearing 
in the testing functional, as well as the use of Schwarz's inequality for 
functions, to eliminate the spatial dependence from the testing functional. 

The difficulty in using this method is that the definition of integral 
coordinates is not always possible. Moreover, the method generally yields 
conservative results. The assumed modes method discretizes the system 
by representing the continuous displacements by finite series of space- 
dependent admissible functions multiplied by time-dependent generalized 
coordinates. Integration over the elastic domains eliminates the spatial 
dependence, so that the testing functional reduces also in this case to 
a testing function. The main criticisms of the method are the truncation 
effect, which generally leads to more conservative stability criteria, 
and the amount of labor involved in deriving the criteria. 

Another aspect of the stability analysis is the definition of equili- 
brium. In certain cases, the equilibrium is one in which not all the 
coordinates are zero (see Ref. 4). In such cases, the equilibrium is 
referred to as nontrivial, and it is necessary first to solve for the non- 
trivial equilibrium and then to expand the testing function about this 
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equilibrium. Note that a typical example of nontrivial equilibrium for 
flexible spacecraft is that in which the flexible parts are deformed 

f 

under centrifugal force. 

This report presents a stability investigation of the GEOS satellite 
by the assumed modes method. In considering the stability of small motions 
about nontrivial equilibrium, it is shown later that if the analysis per- 
formed by ignoring the motion of the mass center indicates stability, then 
the system remains stable if the motion of the mass center is included. 

Derivation of the Testing Functional 

Let us define a set of body axes xyz as the principal axes of the body 
in nominal undeformed state. We shall refer to these axes as a global 
system . In addition, let us define sets of axes x^z.. (i = 1,2,..., 8) such 
that x^ is directed along the length of the elastic members in undeformed 
state and y^ and are perpendicular to x. . The set of axes x^z. will 
be referred to as a local system . The motion of the spacecraft can be des- 
cribed by the rotational coordinates e.(t) (j = 1,2,3) of the global system 

J 

xyz and by the elastic displacements v^(x.|,t) and w^(x^,t) {i = 1,2,..., 8) 
relative to the local system x.y^z. . In general, the displacements of the 
elastic members cause the mass center of the spacecraft to move relative 
to its nominal position, where the latter is identified as the origin of 
xyz. It is shown in Ref. 3, however, that this shift in the position of 
the mass center can be ignored without affecting adversely the stability 
criteria. Moreover, assuming that the mass center of the spacecraft moves 
in a known orbit in space, the kinetic energy of rotation about the mass 
center can be written in the matrix form 
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(1) 
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IT T 1 u • T • 

w vJr,nr 1 . r . « - {u!} {u!}dm. 

m. r i i 


T = ^{a)> , [J]{to} + {to} { K} + J b 2 j, 

where {w} is the column matrix of the angular velocity components and [J] is 
the inertia matrix of the deformed body. Moreover, {K} is the angular mo- 
mentum matrix due to elastic velocities alone and {uj} is the matrix of 
the elastic velocities relative to x^y^z^. 

The potential energy is entirely due to elastic deformations and can 
be written in the form 
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tL L i=l 'o 


r*i 

EI. 

1 

r 2 2 

3 V. ? 3 W. ? 

<V> +< .Vi 

■'0 

3X i SX i 


dx. 


1 E P 

2 i=l Jo X1 


3V_. ? 3w i ? 

(«4 2 * 


dx, 


(2) 


where El. and P . (i = 1,2,..., 8) are bending stiffnesses and axial forces, 

1 XT 

respectively. The functions v.(x^,t) and w^(x^,t) are subject to given 
boundary conditions. Note that in the case of the members 3 and 4 the 
bending stiffness is zero (or nearly zero) and one of the boundary con- 
ditions at x^ = £.j depends on the tip mass m. . 

Because this is a natural system, the Hamiltonian is simply 


H = T + V a (3) 

But the spacecraft is torque-free, so that the angular momentum about the 
mass center must be conserved. It is shown in Ref. 1 that the conservation 
of the angular momentum can be expressed in the matrix form 

[0]{to> + {K} =’ { 0} (4) 

where {$} is the matrix of the conserved angular momentum. Introducing 
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Eq. (4) into Eq. (3), in conjunction with Eq. (1), the Hamiltonian reduces 
to 


where 

*1 8 r . t • 

= A s {u!} {u i }dm. 

c i=1 Jin. 1 1 1 

is a quadratic function of the elastic velocities and 

T o = 


(5) 

( 6 ) 

(7) 


depends on the generalized coordinates alone. Moreover, only two of the 
angular coordinates are present in T Q . This can be easily explained by 

means of the following argument. Assuming that initially the direction of 
the angular momentum vector coincides with the inertial axis Z and that its 
magnitude is 6, then after some perturbation the angular momentum vector can 
be written in the matrix form {$} * where is the column matrix 

of the direction cosines between axis Z and axes xyz. These direction 
cosines can be expressed in terms of only two angular coordinates. 

For the system to be stable in the neighborhood of the equilibrium, 
it is necessary that the Hamiltonian be positive definite (see Ref. 1). But 
T 2 is positive definite by definition, so that H is positive definite if 
the functional 


K 


= T + 
0 


EL 


( 8 ) 


is positive definite. The testing of k for positive definiteness is hindered 
by the fact that k is a functional and not a function, as it involves the 
continuous variables v- and w. in integral form. We shall circumvent this 
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difficulty by using the assumed modes method. 

The term T Q in Eq. (8) involves the motion x c , y c , z c of the mass center. 
Generally, these terms complicate the stability analysis enormously without 
affecting materially the stability statement. Indeed, it is not difficult 
to show that if the system is judged as being stable in the sense of Liapunov 
on the basis of an analysis that ignores x c , y Q , and z Q , then the same con- 
clusion is valid for the actual motion. To this end, let us write 


CO] = [J] u - [J] c 


(9) 


where [J] u is the inertia matrix obtained by ignoring x c , y c , and z Q and 


[0] c = 


• 2 2 
+ z c 

"Vc 

-Vc 


"Vc 
2 2 
4 + z c 

"Vc 


-x„z_. 
c c 

-y z 

J c c 

2 2 

x"; + y~; 
c c 


(10) 


Whereas the matrices [0] and [J] are positive definite, the matrix [J] c 
is only positive. From Eq. (9), it follows that for any arbitrary vector 
{a} the quadratic forms associated with [J] and [J] u satisfy the inequality 


From Appendix B, however, we conclude that 


(ID 


{e} T [or 1 m im T [J],7 1 (e} 


02) 


Next, let us introduce the functional 

■c, = ^B} T tJ]^{6} + V a 03) 

By virtue of inequality (12), we conclude that 

k > k ■ (14) 
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so that, if k-| is positive definite the system is asymptotically stable. 

The preceding statement is true irrespective of the magnitude of x c » y Q 
and z , although when they are large the stability criteria derived by 

w 

using k-j insteady of k can be very restrictive. In most practical cases, how- 
ever, x , y , and z are one order of magnitude smaller than the elastic 
c c c 

displacements themselves, in which case appreciable simplification of the 
stability analysis is achieved by ignoring them, without sacrificing accuracy. 


Calculation of Nontrivial (Deformed) Equilibrium 
a . Problem formulation 

The equilibrium state to be considered is that in which the spacecraft 
spins about the symmetry axis with constant angular velocity Q, as shown 
in Fig. 1. In that state the radial booms remain undeformed, but the axial 
booms undergo bending deformations in two perpendicular directions as a 
result of the centrifugal forces. The distributed centrifugal forces are 
equal to the negative of the distributed mass multiplied by the centripetal 
accelerations . Hence, we wish to calculate first the centripetal accelera- 
tions. Considering boom i (i = 5, 6, 7, 8), we can write the position vector 
of any point on the boom in the form h ^ + r^ + u^, where h. is the vector 
from the satellite center to the point of attachment of the boom, r . is the 
vector from the point of attachment to any arbitrary point on the boom, and 
u i is the corresponding displacement vector. Denoting by i., j . , the 
unit vector along the local axes x., y. , z., the position vectors are as 
follows 


h^+r.+u. = (h xi +x i -z c sina i )i 1 + (h^.+v^j^- + (h 2i +w.-z c Cosa.)k. , i=5,6,7,8 (15) 


Recognizing that = ftk, where k = i ^ s i na . + k^cosa^, the centripetal accelera- 
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tions are 


2 

a. = ^ 2 >c[qx ( + r. + u.)] = n {[(h z .j + - z c cosa.|)sina i . 

-(h x . + x i - z c sina i )cosa i ]cosa i i i - (h y1 + v^j-j 

-E(h zi + - z^osa.Jsina^ - (h x< . + - z c s1na i )cosa.]sina 1 .k 1 .} (16) 

Hence, neglecting the relatively small quantities w.j, the centrifugal axial 
forces become 

P xi = - I p i n[h zi sina i - (h xi + x.)cosa i ]cosa.dx i 

= P i fi 2 {^{(h xi +)!. i ) 2 -(h x .+x.) 2 ]cosa 1 -h zi ()i 1 -x i )sina 1 .}cosa 1 

i = 5, 6, 7,8 (17) 

where p. is the constant mass density. On the other hand, the transverse 

distributed forces are 
2 

p . = p*ft (h • + v * ) 

H yi ' yi i' 

p zi = P^ 2 [(h zi +w 1 )sina i - (h xi+ x.)cosa i ]s1na i i = 5, 6, 7,8 (18) 

(19a) 
(19b) 
(20a) 

w i0 (°) = W ; o (°) - °, w ^(£.) = wj''^) = 0, 1=5,6, 7,8 (20b) 


The differential equations and the boundary conditions for the 
equilibrium deformations (x^) and w. Q (x^) are 


d\. . dv. 

El — - -ji- (p . -5-^) *= p , ’ , i « 5,6, 7, 8 
d x 4 dx.. ' xi dx. J K yi ’ 


v i (0) = v! (0) = 0 , v!^(t.) = ^ * ( «- i ) =0, 1 = 5, 6, 7, 8 


EI ^'^7 {P xi^f )=P zi * ’= 5 ’ 6 ’ 7 - 8 

A i 
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Let the solution of Eqs. (19) have the form 


v io (x i } = 


p 

2 

j=l 


a i j ^ 


1 = 5,6, 7,8 


( 21 ) 


where d>..(x.) are the modes of the fixed-base cantilever beam, in which the 
MO v 

first index denotes the beam number and the second the mode number. The 
explicit expression of <J> ^ j ( x ^ ) is (see Ref. 6, Sec. 5-10) 

♦ ij (x i ) = A id [cos 3 ij £ i + cosh e^.A.Msin P^x. - sinhp^x.) 

- (sin + sinh P-jj^-Kcos p^x. - cosh p.jX.)] (22) 


where the amplitudes A.j are such that the functions are orthonormal, 

i.e., they satisfy relations 


f* 1 ' p, * ik (x i )dx i = « jk 


(23) 


where «- k is the Kronecker delta. Inserting Eqs. (21) and (22) into Eqs. (19), 
we conclude that the coefficients a., must satisfy the algebraic equations 

^ J 


P 

s a 
M 


r, ri. d<f>..d<J>.. rZ- 

U ( “ij 8 jk + L p xi d^dJT dx i } = L Pyi^ik dx i 


xi dx. dx- 'i J 


i - 5, 6,7,8; k = 1,2,. ...p (24) 


Similarly, letting the solution of Eqs. (20) be 


w. 

10 


(x.) = 2 b,. ♦• i (x 1 ), i = 5, 6, 7,8 
• ij 'J 1 


(25) 


we arrive at the algebraic equations 


fa b ij (w ij s jk + 


rt, d*..(i*. k rt. 


i — 5,6, 7, 8; k — l,2,...,p 


(26) 
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to be satisfied by the coefficients b^. 


Pro bl em Discretization by the Assumed-Modes Method 


Next let us transform the functional into a function and, to this 
end, let us derive an explicit expression for T Q . First, we recognize that 
the inertia matrix [J3 U can be written in the general form 

[0] u = [J] r + S [i/wpcip W 

where [J] is the inertia matrix of the rigid hub and [J_j] is the inertia 
matrix of the member i in terms of local coordinates. Its elements are 


J il1 = J p i [(h yi + V io + v n )2 + (h zi + W io + w il )2]dx i 
J i22 = J p i [(h xi + x i )2 + (h zi + W io + w il )2]dx i 

J i33 = f p i C(h xi +x i )2 + (h yi +v io + v il )2]dx i 

d il2 = J i 21 = " , p i^ h xi + x i^ h yi + v io + v iP dx i 

Pi (h x1 + x.)(h z . + w. 0 + w 11 )dx 1 

J i23 = J i32 = ' ) p i (h yi + v io + v il )(h zi + w io + M il )dx i 


J il3 J i31 


(28) 


where v^ Q and w^ Q are the equilibrium elastic displacements and v^-j and w.-j 
are small perturbations. Moreover |>.] is. the matrix of direction cosines 
between the local coordinates x ^ y ^ z * and the global coordinates xyz. It 
will prove convenient to separate the various orders of magnitude in [J] u . 
To this end, let us write 


[J] u = [J] 0 + [03! + [J] 2 


(29) 


where 
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[0] n = z U 1 ] l [J i ] 0 [«- i ] 
0 1=0 ‘ u 


(30a) 


in which 


PL - 


[0 ] - 
L o J o 


A 0 0 
0 B 0 
0 0 C 


(30b) 


is the inertia matrix of the rigid hub in which A, B, and C are the principal 
moments of inertia about x, y, and z, respectively. Note that the above 
statement implies that the global axes xyz are principal axes for the space- 
craft. Moreover, 


<-Wo “ J ^ C(h yi + v io ’ 2 + (l Vi + w io^ dx i 
< J i 22 ^o = J p i C(h xi + x i )2 + (h zi + w io )2]dx i 
< J i 33 ) o = ( p i t(h xi + x i )2 + (h yi + v io )2]dx i 


(30c) 


^ J il2^o = ^ J i2po 


p i (h xi + x i )(h yi + v io )dx i 


< J 113>o = (J i31 } o = ' J p i (h xi + X i )(h zi + “io )dx i 

(J i23 } o = (J i32 } o = ' | p i (h yi. + v io )(h zi + w io )dx i 

where v. and w. are given by Eqs. (21) and (25). Recalling Eqs. (23), we 
TO 10 

can write 


p i v io (x i )dx i = a ij 


p i*i.i (x i )dx i 


p i w io (x i )dx i = * b ij J p i*ij (x i )dx i 

J * 

p i v io (x i )dx i = j, J p i a ij a ik»ij (x iKk (x i> dx i = a ij 


(31) 
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( Pi w io (x i )dx i = b ij 

P 

E 

P 


} p i x i v io (x i )dx i = a iij p 1 x i + ij (x i )dx i 

J Pi x i w io (x i )dx i = ^ b ij J p iVij (x i )dX i 

( p i v io (x i )w io (x i )dx i = J, J, 1 p i a ij b ik^ij*ik dx i 


(31 cont'd.) 


j=l k=l a iJ bik6ik = j ai j blj 


Next, let us write 

[J], = E [^W 
and introduce the generalized coordinates 
6 j = q^(t) , j = 1*2 

p+2 

v n = .H, <i>ij( x 1 ) c, j( t ) 

~ _ (2i-l )p + 2 

v., = I ♦ 1i (x 1 )q i (t) 

11 j=2(i-l_)p+3_ 1J 1 J 

(2n-T)p+2~ 

v, = E W x J q i (t) 

n1 j=2(n-l)p+3 nJ n J 


(32) 


2p+2 


w 


11 


= . * ^i.i (x i )q j (t) 


j=p+3 


21 p+2 

W S1 = t 'l> ii (x.)q i (t) (33) 

11 j=(2i-1)p+3_ 1J 1 J 

2np+2 

W , = E * ,(x )q.(t) 

n1 j=(2n-1 )p+3 nj n J 


Then, we have 
(2i-l ) p+2 

[J.]i = E q,(t) 
1 1 j=2(i-1 )p+3 J 


2 f Pi (h yi +v io )$ ij (x i )dx i ‘I p i (h x1 +x i )<|, ij (x i )dx i 


Pi(hxi +x i)*ij( x i )dx i ° 

V 


2 Jpi (h yi +w io ) *ij (x i )dX i 


2ip+2 

E qj(t) 

j=(2i-l )p+3 J 


0 -Jpi( h yi +W io U ij (x i )dX i 2 } p i (h yi +v io^ij (x i )dx i. 

Pi( h zi +w io)*ij dx i 0 “j p i ( h xi +x i j dx i 

0 2 fp i ( h zi +w 1o )'l' i j( x i ) dx i -| p i (h yi +v. 0 )< | - ij ( x i ) dx i 

- P i (h xi +x i )'l- ij (x i )dx i “jpi ( h yi +v 1 0 j4, ia (x i )dx i 0 

( 34 ) 


so that 
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n T n ( (2i-l)p+2 T v 

[J], - I - .2 , I q i (t)[r.] T [J i ]^ i U i ] 

1 i=l 1 111 1=1 < j=2(i-l )p+3 3 1 i U i 


2ip+2 T 

j=(2i-1 )p+3 J 1 1 1J 1 


( 35 ) 


where [0 - ]Y - and are the corresponding matrices in Eq. (34). Simi- 

1 l J i • J 

larly, we can write 


[j] 2 = I Ca i ] T CJ i ] 2 C^ i ] 


where 


(21-1 )p+2 (2i-1)p+2 

[J,], = E E q l -(t)q k (t) 

1 L j=2(i-1 )p+3 k=2(i-1)p+3 0 K 


Jpi^j( x l)*ik( x T )dx i 0 


(36) 

0 

0 


21 p+2 2ip+2 

+ E E q,(t)q.(t) 

j=(2i-l )p+3 k=(2i-l )p+3 J 


(21-1 ) p+2 21 p+2 

+ E E q,(t)q k (t) 

j=2(i-l)p+3 k=( 21 -1 ) p+3 3 


0 Jp 1 * 1J (x i )* 1k (x i )dx 1 ] 

J p 1+1J (x 1 , *1k {x 1 )dx 1 0 ° 

fp1+1j( x l)+1k (x 1 )dx 1 ° 


,p .A. .(x . (x*)dx. 

J i J i ik v i 


)pi*ij (x i^ik (x i )dx i 0 


(37) 


If we choose the functions ^(x^), ^(x-j) an( ^ ^ij^ x i^ i k ^ x i ^ suc ^ 

f P 1 * 1J (x 1 )* 1k (* i )‘«x i = « jk 

I <*1*1j< x 1>+1k (x 1 )dx 1 = s jk {38) 

I P1*u< x 1>*1k {x 1 )dx 1 = 5 1k 
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then 


(21-1 )p+2 

[JJ? = S q?(t) 
j=2(i-l)p+3 3 


{2i -1 )p+2 

+ . 2 q,-(t)q i+ (t) 

j=2(i-1 )p+3 3 J+p 


1 

0 

0 

2ip+2 , 

1 

0 

0 

0 

0 

0 

+ . , * x qf(t) 

j=(2i-1 )p+3 3 

0 

1 

0 

0 

0 

1 

0 

0 

0 


0 0 0 
0 0 1 
0 1 0 


so that 


n T n ( (21-1 )p+2 0 T 

[0] 2 = .1 - S T. q2(t) U.] T [0.]^.] 

i = l 1 w 1 i=l f j=2(i-l )p+3 J 1 w 1 


2ip+2 
+ 2 
j=(21-l )p+3 


(39) 


qji.ivpowpr j 


(40) 


Note that [J] can contain static elastic displacements caused by centri- 
fugal forces resulting from steady spin 9 whereas [J]^ and [J] 2 contain oscil- 
lations about the deformed equilibrium. To evaluate [0]^ and [J ] 2 , we need 
the matrices [a^]. From Fig. 1, we conclude that the matrices of the di- 
rection cosines are as follows: 


- 


[*■3] 


[* 5 ] = 


[t 7 ] = 



1 

0 

0 



-1 

0 

' 0 




0 

1 

0 

5 

U 2 ] 

0 

-1 

0 




0 

0 

1 



0 

0 

1 




0 

1 

0 



"0 

-1 

0 




-1 

0 

0 

> 

11 

1 — 1 
1 — 1 

l 

0 

0 




0 

0 

1 



0 

0 

1 












— 




0 

cos 

a 5 

sin 


"0 

-COS cv 

sin a g 


-1 


0 

0 

’ [^5] 

1 


0 


0 


0 

r— 

-sin 

8r 

5 

COS a r 
5 


0 

sin ag 

COS ci^ 


cos 

a 7 

0 

sin 


-cos 

a 8 

0 

sin ag“ 



0 

1 

0 

■ C* 8 ] = 

0 


-1 

0 

l-sin 

a 7 

0 

COS Cty 


sin 

a 8 

0 

cos ag 


(41) 


and note that = a, and a, = a„. 
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Assuming that [0]-j and [J] 2 are small compared to [J] Q , we can write 


the following approximation for the inverse of 


[0]' 1 = [K] = [K] o + [K] 1 + [K] 2 


(42) 


where the subscripts 0, 1 and 2 once again identify the order of magnitude 
of the quantities involved, in which 


[K] 0 = [J]; 1 

[K] 1 - -CJ]" 1 

[k] 2 = -[J]; 1 tJ3 2 [J]; 1 + CJ3^ 1 [J]iC03“' , [j] 1 C0 ]- 1 


(43) 


From Eqs. (35) and (43), we conclude that 

n ( (2i-1)p+2 T 

[K] 1 = - [K] [J] [K] = - Z E q (t) [K] [t ] T [J ]y Jt.lCKl 

1 oio i=1 ( j=2(i-1 )p+3 J 01 1 13 1 1 

2ip+2 


+ j=(2i-l)p+3 


q j (t)[K] o [t ;i ] T [0 i ]” j [£ 1 ][K] o j (44) 


Next, let 




(45) 


Then 


CK], = 


n ( (21-1 )p+2 2i'p+2 

M , Z q^t) [A.],.[K] + l q..(t)[B,], .[K] 
i=l 'j=2(i-l )p+3 J 1 |J 0 j=(2i-l )p+3 J 1 lj 1 


(46) 


Moreover, using Eqs. (40) and (43), we can write 

[K] 2 = - [K] o [0] 2 [K] o + CK] o [JD 1 [K] o [0] 1 CK] o 

n ( (2i-l )p+2 2 Tv 

- - z Z q^(t)[K] u ] T [d ]"[£ ][K] 

i=1 * j=2(1-l )p+3 J 01 1 1 1 0 


I 
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2ip+2 0 T 

+ I q.(t)[K] 0 [it ] T [J ]«U ][K1 

j=(2i-1 )p+3 3 01 1 1 1 0 

(2i-l)p+2 T 

+ I q.(t) qi ..Jt)[K] Tz ] T [d ]^ w [. ][K] 

j=2(i-l )p+3 J J+p 0 1 1 2 1 ' 


!’( 


(2i-l )p+2 2ip+2 \) d 

2 + j={2i .^ p+3 CK] < 


=2 ( i -1 )p+3 


( 47 ) 


so that, letting 

ca 1 ] 2 = [Ky^fcj^^u.] 

^] 2 = [K] 0 Ct i ] T [O i ]«[t i ] 
[C.] 2 = [K] 0 [A i ] T tJ i 32 W C^ i ] 


(48) 


we obtain 


n ( (2i-l )p+2 , 2i p+2 , 

DO, = ■ I l q-(t) [A.]o[K] + i q?(t) [B.]„[KL 

4 i=1 (j=2(i-i p+3 J Wo j=(2i-i )p+3 J 12° 


(2i-l)p+2 

+ , ? q i ( t )q 1+D (t) 

j=2(i-l )p+3 3 J+p 120 


n n i (2i-l)p+2 (2j-l)p+2 

+ t t S S q( ,(t)q m (t)EA.] 1 [A.], [K] 

i=l j=l ( n=2(i-l)p+3 ra=2(j-l )p+3 *• m 1 '*■ J lra ' 


21 p+2 

+ z 
£=(2i-1 )p+3 m=(2j-l)p+3 


2 j p+2 

t q,(t)q m (t) [B.] u [B.] lm [K], 


(2i-l)p+2 2jp+2 

+ 2 I q„(t)q m (t) [A-], [B .], [K] 

t=2(i-l )p+3 m=(2j-l)p+3 1 m i U j lm o 


21 p+2 (2j-l )p+2 

+ s t q,(t)q (t) [B.], [A,], [K] 

£.={21 -1 ) p+3 m=2(j-l)p+3 1 m i U J lm o 


(49) 


Introducing the notation 



Eqs. (46) and (49) can be rewritten in the form 

n ( (2i-l)p+2 2ip+2 

[K]] = - £ j £ q,(t) [R.]- + £ q.(t) [SJ,. 

. 1=1 <j=2(i-1 )p+3 J 1 J j=(2i-l)p+3 J 1 J 


( 51 ) 


and 


n r (2i-l)p+2 9 2ip+2 9 

[K] ? = - M 1 [RJ ? + s qft) [S.] ? 

i=l I j=2(i-l)p+3 3 1 c j=(2i-l)p+3 3 1 1 


(21-1 )p+2 

+ £ q.(t)q. (t) [T.-L 

j=2(i-l )p+3 3 3 p 12 


+ small terms 


(52) 


Next, we wish to write the expression for k . Assuming that the orien- 
tation of the global system xyz is obtained from the inertial space XVZ by 
the rotations about z, 6-| about x, and about y, then the direction 
cosines between axes xyz and Z are as follows: = -cose^sin02, n ^ ” 

sin0-|, = cos0 1 cose 2 . For small angles e-| and the column matrix 

iij) can be approximated by 


U z > - (*7 } o + ^z 1 ! + {5 -z } 2 


(53a) 


where, recalling that = q-j and we have 


<*Z } o = < 0 



U Z } 2 




(53b) 


Recognizing that s = Csi, the functional k, can be approximated by 
K 1E = 2 c2 n 2 {{£, z >J[K] o {K. z > 1 + 

+ 2(«. z >{[K] 1 {t z } 0 + {t z } J[K] 2 { «. z } 0 ) + % (54) 


in which 
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tVl^Wl ^ K ldo q 2 + ^ K 22^o q ? + 2 ^ K 12^0 q l q 2 

^A^Wo ’ -< K 33>o (q l + q 2> 

2{S 'z } l[K]l {!! 'z } 0 = -2 q 2 ( ^ 3) 1 + 2q 1 (K 23 ) 1 

n f (2i-l )p+2 

= ' 2 i = 1 lj=2 {i -l)p+3 q j (t) [ " (R il3 ) lj q 2 + (R i23 ) 1j q l ] (55) 

2ip+2 \ 

+ j=(2i-l)p+3 qj(t) C ' (S il 3) lj q 2 + (S i2 3) lj q l ] ) 

T n (2i-1)p+2 ? 

‘ {£ 7 ) n [K] ? {A } = - 2 2 q^(t)(R.-o)o 

2 0 2 z 0 i=l ^j=2(i -1 )p+3 J 13 2 


. 2ip+2 


( 2 i-l )p +2 


2) p • (t) {S .- 3)2 * £ P-; (t)p -?+ n (t)(Tj 33) 0 

j=(2i-l)p+3 3 133 2 ,i=2(i -1 )d+ 3 3 3+p 133 2 


j=2(i-l)p+3 


But s by virtue of the fact that the functions and $ . . satisfy corresponding 

* J ^ J 

eigenvalue problems, the elastic potential energy satisfies the inequality 

(56) 


V EL > " 


1 " ( < 2i ‘ 

2 i=l V j= 2 (i- 


-l)p +2 « ? 2 ip +2 

2 At .q* + 2 A 

l)p+3 J k=(2i-l )p+3 


2 „ 2 ) 

ik q k/ 


[H] = 


where A., and A*, are the natural frepuencies associated 'with the modes <j>.. 

1 j 1 K I J 

and Replacing in Ep. (54) by the expression on the right side of 

inequality (56), the system can be regarded as asymptotically stable if 

(57) 

is positive definite where [H] is the Hessian matrix given by 

( l< '22 )o - ( ^33 ) o ^12*0 "( R i23hj *” ^ S i2 3 Hk 

^ R i13hj •" ^ S i 1 3^ i k 


K 2E = ? C 2 ! 2 2 {q} T [H]{q} 


'( R n 2 3 h«. 


“I 

n23 / lm 


^ R nl 3^ 1 s, ( S nl3hm 




symmetric 


(A*k> 




^*m> 


n n 


( 58 ) 



where 


(A*j )' 


(A tk ) 


,2 

2 2 ~ ^ ■* o o ) 1 


C^‘ 


‘133 ^ 2 


,2 

-IK- . ($ ) 

cV 1332 


/ j = 2(i-l)p+3, 
l k = (2i-l)p+3, 


(2i-l )p+2 
2i p+2 




(A* )' 
nm 


n& 

C 3 R 2 


nm 


C 2 ft 2 


^ R n33^2 

^ S n33^2 


£ = 2(n-l)p+3, 
m = (2n-l)p+3, 


(2n-l )p+2 
2in+2 


(59a) 


(59b) 


The function is positive definite if the matrix [H], which in turn requires 
that all the eigenvalues of [H] be positive, A computer program has been 
written for the calculation of [H] and for the evaluation of its eigenvalues. 
The program is described in the next section. 


Description of the Computer Program 

The computer program follows in detail the equations derived for Hessian 

matrix. Some explanations of all its subroutines are given below: 

1. The standard Gauss-Jordan is used for the inversion of matrices. The 
corresponding subroutine is called MI NV . 

2. Subroutine HSBG reduces an n by n real matrix A to an upper almost 
triangular form by a similarity transformation. Each row is reduced in 
turn, starting from the last one, by applying right elimination matrix, 
and similarity is achieved by also applying the .left inverse transforma- 
tion. Thus the eigenvalues of A are preserved. Similarity transfor- 
mations are using elementary elimination matrices with partial pivoting. 

3. Subroutine ATEIG computes the eigenvalues of a real upper almost triangu- 
lar matrix (Hessenberg form) using the double QR iteration of Francis 
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(Ref. 7). If all the eigenvalues of the matrix are positive, then the 
matrix is positive definite. 

4. Subroutine GMPRD is used to multiply two general matrices to form a 
resultant general matrix. 

5. Subroutine GTPRD is used to premultiply a general matrix by the transpose 
of another general matrix. The transpose of A is not actually calcu- 
lated. Instead, elements of matrix A are taken columnwise rather than 
rowwise for postmultiplication by matrix B. 

6. DRTMI determines a root of the general nonlinear equation f(x) = 0 in 
the range of x from x^. up to x^ (x^., x^.. given by input) by means of 
Mueller's iteration scheme of successive bisection and inverse parabolic 
.interpolation. The procedure assumes f(x .)*f(x„.) £ 0. Convergence 

is quadratic if the derivative of f(x) at root is not equal to zero. 

All the subroutines described above could be found in the System/360 
Scientific Subroutine Package. 

The function SIMPS is used to evaluate numerically all the integrals 
by n repeated applications of Simpson's rule, where n is given by the NASR 
variable in the program. Because all the chosen admissible functions 
involve only well-behaved curves, use of Simpson's rule for all the inte- 
grations is justified. 

Function FCT contains the equation cos $x cosh $x + 1 which is the 
characteristic equation of the nonrotating cantilever beam. 

Function THI is used to calculate the value of an admissible function 
for either beam or cable with given amplitude, frequency and argueinent. 

Function DTHI is the derivative of function THI and function TTHI is 
the integral of function THI. 

Total moment of inertia [J] and the moment of inertia [ J ^ and 
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for each beam and each cable are calculated numerically by subroutine 
NFREQ. One identification variable (ID) indicates that the input data belongs 
to either a cable or a beam. A corresponding procedure is used to determine 
the coefficients of the admissible functions associated with the nontrivial 
equilibrium. 

Function DIJ and its entry functions provide some products of the 
independent variable x, the admissible functions and their derivatives, 
which are all involved in the centrifugal terms appearing in the differen- 
tial equation of nontrivial equilibrium. 

Finally, all the elements of Hessian matrix are obtained numerically 
in the main program and the property of the matrix is tested by solving all 
the eigenvalues of the matrix using the QR iteration method. 

Numerical Results 

The preceding computer program has been used to test the stability of 
the GEOS satellite (the "simple model"). The numerical data (per letter 
of Dr. Peter Kulla dated 9 May 1974) is as follows: 

*1 = s 2.66m , h xl = h x2 = 0.73m , h ^ = h^ 2 = 0 , h zl = h z2 = -0.5m 

_] 

a-j - a 2 - 0 , p.j = p 2 = 1.127 kgm ; w-j = ta 2 = 2 Hz ; 

&3 ” ^4 _ 20m , h x 2 ~ h x ^ - 0.73m , h^ 2 = h ^ = 0 , = h z ^ ~ 0.15m 

a 3 ~ a 4 = 0 » P3 = p 4 = 0.03 kgm”'' , - m^ = 0.1kg; 

£5 - ~ - ^5 ~ sin27° + 0.5 cos27° m , h^g = h^^ ~ -0.42m 

h z5 = h z6 = 0,8 cos27 ° “ sin27° , c< 5 = - 27° > P5 = P5 = 0.733kgrrf^ 
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ti)g “ Wg “ 3Hz , 


A 7 “ ^8 = 1,5m * h x7 = h x8 = 0,8 sin45° + 0.5 cos 45°m , h y7 = = -0.42m 


h z7 = h z g = 0.8 cos45° - 0.5 sin45°m, ay = ag = 45° , py = p g = 1.132kgm 


-1 


)y “ Ug 5Hz j 

t*- 


A - r xx + 2 


3p 3^ h z3 + ^ h x3 + x 3^^ dx 3 + 2m 3^ h z3 + ^ h x3 + * 3 ^ = 354.13806kgm‘ 


B = I = 125kgm' 


C = I +2 
zz 


ft - 1 rad s 


Jo 

-1 


^ p 3^x3 + x 3^ x 3 + ^ m 3^x3 + ^3^ = 397.10606kgm^ 


First, the nontrivial equilibrium configuration was evaluating by using 
two terms in series (21) and (25). The results are as follows: 


a., = -0.13710 x 10' 2 mkg 1/2 

a. - = -0.19418 x 10' 4 mkg 1/2 

b il = -0.36107 x 10' 2 mkg 1/2 

b. 2 = -0.22324 x 10‘ 4 mkg 1/2 

a., = -0.43397 x 10" 3 mkg 1/2 
a i2 = -0.61305 x 10' 5 mkg 1/2 

b^ = -0.92787 x 10‘ 3 mkg 1/2 
b- 2 = -0.74565 x 10" 5 mkg 1/2 


v i - (£..) = -0.18229 x 10' 2 m 

w. (£.) = -0.48396 x 10" 2 m 
i = 5.6 

v io (£ i ) = -0.65666 x 10" 3 tn 

w-qC^i ) = -0.14127 x 10' 2 m 
i = 7.8 


Using the above results, and using two terms in the series (33), a 
34 x 34 Hessian matrix was obtained. The matrix failed the test of 
positive definiteness, a fact that can be traced to the cables. 
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Conclusions 


The stability problem associated with a spin-stabilized satellite 
similar in configuration to the GEOS satellite has been formulated and 
programmed for digital computation. The formulation is capable of 
accommodating satellites with a somewhat different configuration than the 
GEOS, in the sense that the number of elastic members and their orientation 
relative to the spacecraft is arbitrary. 

For the given configuration, one eigenvalue was found to be negative, 
so that on the basis of the Liapunov direct method the spacecraft cannot 
be judged as being stable. By inference, the system can be regarded as 
being unstable. This lack of stability can be traced to the fact that the 
lowest natural frequency of in-plane vibration of the cables is close to 
zero (see Appendix A). This is based on the assumption that the 
cables do not possess bending stiffness. In view of the negative stability 
statement obtained, a study of the effect of small cable bending stiff- 
ness on the spacecraft stability appears warranted. However, even if an 
analysis including small cable bending stiffness indicates stability, the 
first natural frequency of oscillation of the spacecraft is likely to be 
very low, so that the cables can represent a potential problem area. 
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Appendix A - Eigenvalue Problems for Rotating Elastic Members 


a. Rotating Beam 

Let us consider a rotating beam in transverse vibration, as shown in 
Fig. 2. The eigenvalue problem is defined by the differential equation 
(see Ref. 6) 

EI - \ pn 2 ^ {[(h+Jt) 2 - (h+x) 2 ] ^ } = A 2 p$ , 0 < x < £ (A.l ) 

where (x) is subject to the boundary conditions 

4> = 0 and = 0 at x » 0 (A. 2) 

2 3 

— i = 0 and = 0 at x = £ (A. 3) 

dx^ dx J 

There is no closed-form solution of the eigenvalue problem (A.l) - (A. 3). 
Hence, we wish to obtain an approximate solution. To this end, we use the 
Rayleigh-Ritz method and assume a solution in the form of the series 

n 

<f>(x) = z a u (x) (A. 4) 

r=l r r 

where u r (x) are comparison functions, namely, functions satisfying all the 
boundary conditions of the problem but not the differential equation (otherwise 
they would be eigenfunctions). We choose as comparison functions for the ro- 
tating bar the eigenfunctions of the nonrotating cantilever beam, obtained 
by setting 0 = 0 in Eq. (A.l). These functions are (see Ref. 6) 

u r (x) = A r [(sin$ r L - sinh$ r L) (sine r x - sinh 3 r x) 

+ (cos3 r L + cosh$ r L)(cos6 r x - cosh0 r x)] , r = l,2,...,n (A. 5) 

where the coefficients A^ are arbitrary $ r are the eigenvalues of the 
problem; they satisfy the characteristic equation 


A-l 



cos$ r L coshe r L +1=0 


(A. 6} 


The functions u^(x) are orthogonal. Moreover, it will prove convenient to 
remove the arbitrariness from u r (x) (r = l,2,...,n) and determine the 
coefficients A r uniquely by normalizing the functions u r (x) so as to satisfy 

f £ 

pu r (x)u s (x)dx = s rs , r,s = 1,2 n (A. 7) 

where <$ rs is the Kronecker delta. 

It can be shown that the Rayleigh-Ritz method, in conjunction with the 
normalized comparison functions u r (x), lead to the special eigenvalue problem 

[k]{a> = A 2 (a> (a. 8) 


where the matrix [k] is real and symmetric; its elements have the values 


2 ? 
d u d u 


r* a u a u i o o d 

k - El — * dx + i p ST [ ( h+Jt ) ^ - (h+xr] 3 

rS J 0 dx^ dx^ L In d 


du^ du n 
r s_ 

dx dx 


dx 


2 1 2 
= w 6 + 4 pft 

r rs 2 K 


f £ ? du du 

[(>m) - (h+x) ] gJlg-idx , r,s = 1,2 n (A. 9) 


in which y are the natural frequencies of the nonrotating beam. The solution 
of the eigenvalue problem (A. 8) and (A. 9) consists of the eigenvalues A?, which 
are the squares of the estimated natural frequencies of the rotating beam, and 
the eigenvectors {a^h (i = l,2,...,n). It follows that the estimated 
eigenfunctions are 


♦|(x) 


n 

s 

r=l 



(A. 10) 


Eigenvalue problems of the type (A. 8) and (A. 9) must be solved for 
radial members such as 1 and 2. For members 5, 6, 7, and 8 the eigenvalue prob- 
lem must be modified to account for the inclination of the bar and the resulting 
transverse loads. 



b. Rotating cable with tip mass 

The eigenvalue problem for a rotating cable is similar to that of the 
rotating beam shown in Fig. 2, except that the bending stiffness is equal to 
zero. In addition, we are interested in the case in which the cable has a 
tip mass m. The corresponding eigenvalue problem is defined by the differential 
equation 

- fi 2 pEOvH) 2 - (h+x) 2 ] + m(h+£) j = A 2 p* , 0 < x < p (A. 11) 

where <j>(x) is subject to the boundary conditions 


* = 0 at x = 0 (A. 12) 

-m(hH)ft 2 mA 2 4> = 0 at x ~ £ (A. 13) 

The eigenvalue problem (A. 11) - (A. 13) has no closed-form solution either. 
The eigenvalue problem of the rotating string with no tip mass and with h = 0, 
however, is satisfied by the Legendre functions. The Legendre functions of odd 
degree can be used as admissible functions for the eigenvalue problem (A. 11) - 
(A. 13) as they solve a similar problem and satisfy the boundary condition at 
x = 0. Note that admissible functions need satisfy only the geometric boundary 
conditions of the problem. Hence, let us assume an approximate solution in 
the form 

n 

*(x) = i a P„ ,(x) (A. 14) 

r=1 r dr ~' 


where 

Pl(x)=f 

P 3 (X) = 1 [5(f) 3 - 3 f] 


P 5 (x) = 1 [63(f) 5 -.70(f) 3 + 15 f] 


fl-3 


(A. 15) 


(A. 16) 


are known as Legendre polynomials. They possess the orthogonal ity property 
0 

and they satisfy the relation 
rt 


ft 

j *V X P k^ dx “ 0 ’ j » k = 1 , 2 ,... 


0 


P j (x)dx = 2jTT ’ 


(A. 17) 


The Rayleigh-Ritz method leads to the eigenvalue problem 

[k]{a} = A 2 [m]{a} (A. 18) 

where the matrices [k] and [m] are real and symmetric. Their elements are 
ri 
0 


k rs = 


{^p[(h+£) 2 - (h+x) 2 ] + m(hn)}P' r _ 1 (x)P^ s _ 1 (x)dx 

r,s = 1 >2,. . .,n 


(A. 19) 


and 

m rs 


= p j P 2r-l (x)P 2s-l (x)dx + mP 2r-l (l)P 2s-l W 


” 2(2r-l )+l 6 rs + mP 2r-l ^ P 2s-1 ^ (A. 20) 

The solution of the eigenvalue problem (A. 18) - (A. 20) consists of the 
2 

eigenvalues , which are the squares of the estimated natural frequencies of 
the rotating string with a tip mass, and the eigenvectors {a^} (i = 1 , 2 ,... 
It follows that the estimated eigenfunctions are 

*,-(*) = * a r l)p 2r-l (x) (A. 21) 

The eigenvalue problem (A. 18) - (A. 20) must be solved for member 3, 
which yields automatically the solution also for member 4. 
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Appendix B - Theorems on Inequalities for Quadratic Forms 


Theorem. Given two matrices A and B which are symmetric and positive 
definite over real number field R. If x^Ax >_ x^Bx for any vector x over 
R, then x T A _1 x <_ x T B"^x. 

Proof: Because A is symmetric, there is an orthonormal matrix U such that 

A* 1 ' 2 V UA-W (B.l 

where A is a diagonal matrix with its elements equal to the eigenvalue of 
the matrix A. The effect of the operation 


C = a“ 1/,2 ba“ 1/2 

is to transform the symmetric and positive definite matrix B 
C which is also symmetric and positive definite, namely, 


(B. 2) 

into a matrix 


C T = (A-^BA- 1 / 2 ) 7 = A^Va' 1 ' 2 = A^BA' 1 ' 2 = C 


(B.3) 


Similarly, there exists an orthonormal matrix V such that 

V T CV = vV 1/2 BA' 1/2 V = y (34) 

where n is a diagonal matrix. 

Introducing the linear transformation 

P = vWx ' _ . - (B.5) 

T T 

into the inequality x Ax > x Bx, we obtain 


p T V T A- 1 / 2 AA-^ 2 Vp > pVa^^BA ' 1 / 2 


Vp 


which reduces to 


T, T 
P IP i P up 


(B.6) 

(B.7) 


B-l 



where I is the identity matrix. Because A and B are positive definite, all 
the elements of the diagonal matrix p are positive. It follows from 


inequality (8.7) that 

pVVipVV (B.8) 

Moreover, recalling inequalities (B.6) and (B.7), it follows that 

pVa^VV^Vp < p T V T A 1/2 B _1 A 1/2 Vp (B.9) 

Next, let 

y = A 1/7 Vp = A 1/2 VV T A 1/2 x = Ax (B.10) 

so that inequality (B.9) reduces to 

yVy < yVy . . (B.ll) 


Because A is symmetric and positive definite, we can show that A can be 

regarded as a linear transformation mapping the linear space into itself. 

T -1 T -1 

This concludes the proof that x A x <_ x B x. 

Corollary. Given two matrices A and B which are symmetric and positive 

T T 

definite over real number field R. Then x Ax >_ x Bx for any vector x over 
R if and only if every eigenvalue p^ (i = l,... s n) of A~"*B is such that 

1 _> p.j > 0. 

Consider the series 

I + A^B + (A -1 B) 2 + ... + (A” 1 B) m + ... (B. 12) 

1 m 

For convergence it is clearly necessary that lim (A" B) -> 0. This 

, m m 

condition is also sufficient, for if lim (A" B) m -* 0 if follows that |p r | < 1, 

I m^°° t 

and therefore I - A" B does not vanish and (I-A~ B)~ exists. But 



1 _1 

so that postmultiplication of Eq. (B.13) by (I - A B) yields 

I + A _1 B + (A _1 B) 2 + ... + (A -1 B) m = (I - A^B)' 1 - (A' 1 B) m+1 (I - A _1 B) _1 (B.14) 

As m ■> «, Eq. (B.14) reduces to 

I + A _1 B + (A _1 B) 2 + ... + (A’^B) m + ... = (I - A^B)" 1 (B.15) 

But postmulti plication of Eq. (B. 1 5) by A~^ gives 

A -1 + A^BA" 1 + (A -1 B) 2 A _1 + ... + (A"^B) m A”^ + ... = (I - A^B)"^" 1 = 

[A(I - A^B)]" 1 (B.16) 

Hence, 

A' 1 + A^BA' 1 + (A^BjV 1 + ... + (A^B/V 1 + ... = (A - B)" 1 (B.17) 

so that if two symmetric and positive definite matrices A and B satisfy the 
T T 

inequality x Ax > x Bx for any vector x over R, then the series expansion 
(B.17) is valid and convergent. 
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